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1.  Introduction 


The  systems  carrying  electric  charges  are  widespread  in  nature  and  applications. 

Variational  principles  and  methods  play  a  key  role  in  modeling  these  systems,  and 
they  have  been  developed  for  a  couple  of  centuries. 1-5  The  issues  of  stability  and 
instability  of  electrostatic,  magnetostatic,  and  quasi-static  systems  are  extremely 
important  in  multiple  electrostatic  systems  (e.g.,  Eamshaw  theorem  or  the  Rayleigh 
surface  instability)  and  in  multiple  magnetostatic  systems  (e.g.,  the  well-known 
surfacial  patterning  in  ferrofluids  demonstrated  in  Fig.  1). 


Fig.  1  Ferrofluid  on  a  reflective  glass  plate  under  the  influence  of  a  strong/steep  magnetic 
field.  (Reprinted  with  permission  of  Gregory  Maxwell,  Wikipedia: 
https://en.wikipedia.org/_wiki/  File:  Ferrofluid_Magnet_under_glass.jpg#/.) 

Therefore,  it  is  quite  surprising  that  the  key  instruments  in  the  analysis  of  stability 
(e.g.,  second-energy  variation)  are  essentially  underdeveloped.  Also,  various 
variational  approaches,  closely  associated  with  the  first-  and  second-energy 
variation,  require  quite  sophisticated  mathematical  techniques  to  handle  them. 

Variational  approaches  are  of  key  importance  in  computational  implementations  of 
the  models.  It  happens  because  of  the  a  priori  guaranteed  mathematical  self- 
consistency  of  these  methods  and  open  possibilities  of  using  powerful  numerical 
schemes  based  on  minimization  and  optimization.  Nonetheless,  even  a  quick  glance 
at  the  literature  shows  that  the  variational  approach  is  definitely  underdeveloped. 
In  particular,  in  terms  of  calculus  of  variations,  all  of  the  studies  are  based  on  the 
analysis  of  the  first-energy  variations  (i.e.,  analysis  of  ponderomotive  forces  and 
conditions  of  equilibriums.)  The  analysis  of  the  second  variation,  the  cornerstone 
of  the  stability  investigation,  is  not  even  touched  in  those  famous  publications. 
Moreover,  the  very  expressions  of  the  second-energy  variations  still  remain 
unknown  (Grinfeld5-7  and  references  therein  are  very  rare  exceptions  to  this  rule). 

All  of  the  discussed  issues  lead  to  the  conclusion  that  numerical  modeling  remains 
the  only  realistic  chance  of  exploring  various  morphological  patterns  appearing  at 
the  late  stage  of  destabilization. 
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In  the  following,  we  present  an  iteration  scheme  of  possible  numerical  analysis  of 
the  poststability  patterns.  The  approach  is  based  on  the  variational  formulation  of 
the  model  system  under  study.  The  presentation  is  limited  to  discussion  of  the 
classical  Rayleigh  instability  of  a  charged  conducting  liquid  drop.  The  approach 
permits  extension  to  several  other  systems,  including  those  accumulating  and 
carrying  free  electric  charges  not  only  on  the  surface,  but  also  inside  the  bulk,  as 
well  as  the  pyroelectric  crystals  growing  on  rigid  substrates,  patterns  of 
ferromagnetic  liquids,  and  the  like.  To  make  the  presentation  as  simple  as  possible, 
we  limit  ourselves  to  incompressible  substances  and  ignore  polarization  effects. 

2.  Electrocapillary  Instabilities 

Rayleigh  first  predicted8  that  an  isolated  charged  conducting  drop  is 
morphologically  unstable  if  the  electric  charge  is  sufficiently  great. 

The  intuitive  explanation  (not  the  professional  proof)  for  the  Rayleigh  instability  is 
rather  straightforward.  The  electrons  repel  each  other  and  try  to  stay  as  far  as 
possible  from  each  other.  Therefore,  they  create  the  forces  deforming  the  drop.  The 
assumption  of  the  liquid  incompressibility  does  not  allow  the  spherical  drop  to 
change  its  radius.  Therefore,  the  available  shape  modifications  are  not  radial.  One 
of  the  simplest  nonradial  modifications  is  the  spheroidal  perturbation,  as  shown  in 
Fig.  2. 


Fig.  2  Rayleigh  instability  of  a  charged  conducting  incompressible  liquid.  For  a  sufficiently 
high  charge  or  small  surface  energy,  the  spherical  drop  becomes  unstable  with  respect  to 
spheroidal  perturbations. 

From  the  energetic  point  of  view,  we  would  say  that  electrons  try  to  diminish  the 
accumulated  electrostatic  energy.  To  put  it  another  way,  the  electrostatic  energy  of 
the  system  can  be  lowered  by  letting  the  drop  deviate  from  spherical  shape.  Such 
deformation,  however,  is  opposed  by  surface  tension  since  the  minimum  surface 
energy  is  attained  when  the  shape  is  spherical. 

The  overall  instability  criterion,  then,  reflects  the  balance  between  these  forces.  If 
R  is  the  radius  of  the  drop,  Q  is  the  net  electric  charge,  and  a  is  the  surface  energy 
density  of  the  boundary  interface,  then  Rayleigh’s  criterion  for  morphological 
instability  is 
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(1) 


Ra 


Q~ 


16  naR 


>1  , 


which  can  be  rewritten  as 


2 

T 


4  aR 


>1  , 


(2) 


where  r  is  the  surface  charge  density. 

Per  the  Rayleigh  criterion  (Eq.  1),  as  long  as  the  total  charge  Q  or  the  surface  tension 
a  is  large,  the  surface  tension  influence  dominates,  and  the  spherical  shape  appears 
to  be  stable  with  respect  to  small  nonspherical  perturbations  of  the  shape.  But  when 
the  dimensionless  number  Ra  exceeds  the  critical  value  Ra  =  l ,  the  spherical 

shape  becomes  unstable  and  yields  the  spheroidal  shape.  For  example,  a  drop  of 
water  1  mm  in  radius  is  unstable  if  it  is  charged  to  several  volts. 

Rayleigh’s  original  study  was  picked  up  by  many  researchers  including  Tonks,9 
Frenkel,10  Melcher  and  Taylor,11  and  others  (see  extensive  lists  of  relevant 
references  in  Rosensweig12).  Their  research  showed  that,  in  general,  as  in 
Rayleigh’s  original  work,  the  presence  of  sufficiently  large  electrostatic  or 
magnetostatic  fields  leads  to  the  morphological  destabilization  of  plane  interfaces. 
This  class  of  instabilities  plays  a  role  in  the  physics  of  vacuum  discharges  and  in 
the  operation  of  liquid  metallic  emitters  of  charged  particles.  These  phenomena 
continue  to  be  in  a  state  of  active  research.  One  of  the  primary  questions  of  interest 
is  that  of  the  equilibrium  shape  of  destabilized  interfaces.13  From  the  mathematical 
point  of  view,  such  problems  of  equilibrium  are  challenging,  deeply  nonlinear 
problems  with  unknown  boundaries,  reminiscent  of  the  classical  problem  of  the 
equilibrium  shape  of  a  crystal. 

3.  Algorithm  of  Computer  Modeling  of  Equilibrium  Shape  and 
Poststable  Evolution 


3.1  Energy  Minimization  Approach 

The  Rayleigh  problem  can  be  treated  as  the  problem  of  minimization  of  the  total 
energy  W  comprising  2  ingredients,  the  total  electrostatic  energy  and  the  total 

surface  energy  : 


W  =  E  ,  +E 

elec  surj 


(3) 
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The  surface  energy  density  per  unit  area  in  solid  substances  or  in  liquid  crystals  is 
characterized  by  its  dependence  on  the  normal  N  :  <r  ( /V ) .  For  isotropic  liquid  we 
have  just  a  positive  constant  function  cr(jv)  =  a  =  const . 

The  total  electrostatic  energy  depends  upon  the  charge  loading  mechanism.  In 
choosing  the  charge  loading  mechanism,  we  first  supply  a  certain  amount  of  the 
electric  charge  Q  and  then  disconnect  the  battery.  Then  we  arrive  at  the  following 
expression  of  the  total  energy  ingredients: 

Ee,ec-m  I  E,urf^dZ-  <4> 

Space  a 

The  second  integral  in  Eq.  4  is  the  total  area  of  the  conductor’s  boundary  s . 

There  are  actually  2  functional  degrees  of  freedom  in  the  minimization  problem:  1) 
the  unknown  surface  distribution  of  electric  charge  r  and  2)  the  position  of  the 
boundary  h  . 

In  this  minimization  problem,  the  condition  of  mass  conservation  should  be  taken 
into  account: 


J  dZz  N.  =  3Vq  =  const  ,  (5) 

where  n.  is  the  unit  normal  to  the  boundary. 

Also,  given  the  total  electric  charge  Q  of  the  conductor,  we  have  to  respect  the 
constraint: 


\dEr  =  Q  =  const  .  (6) 

It  should  be  taken  into  account  that  the  distribution  of  electric  charge  r  appears  in 
the  total  energy  implicitly.  Indeed,  given  the  probe  surface  charge  distribution  r 
and  the  interface  position,  we  have  to  calculate  the  corresponding  electrostatic 
energy  C(,lec  .  To  do  this,  we  have  to  calculate  the  electrostatic  potential  q> .  In  other 

words,  one  has  to  solve  the  electrostatic  problem  for  that  geometry  and  that  surface 
charge  distribution.  That  goal  can  be  accomplished  by  calculating  the  singular 
integral 


<p(z)  =  J  dZ 


(7) 
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where  z  is  the  radius  vector  of  the  points  on  the  boundary.  The  same  goal  can  be 
achieved  or  by  solving  the  boundary  value  problem  of  electrostatics: 

1)  inside  conductor 


(D  =  (D.  =  const  .  (o) 

T  in 

2)  outside  conductor 

v'vp  =  0 .  (9) 

3)  at  the  interface  liquid-vacuum 

(Pin=(Pout’Vi(PoutNl  =4^-  (10) 

4)  at  infinity 

p->0.  (11) 

Thus,  from  the  mathematical  point  of  view,  we  arrive  at  the  problem  of 
multidimensional  calculus  of  variations  with  the  isoperimetric  and  differential 
constraints,  and  with  the  unknown  boundary. 

3.2  Two-Step  Algorithm  of  the  Energy  Minimization 

Now  we  can  formulate  the  suggested  2-step  numerical  algorithm.  Let  s  ,  r  ,  b  , 
and  w  be  the  current  positing  of  the  boundary  and  the  current  values  of  the  surface 
charge  density,  the  volume  of  the  drop,  and  the  total  energy  of  the  system, 

Hs  *  *  * 

respectively.  We  use  the  notations  s  ,  r  ,  B  ,  and  wun  for  the  same  quantities  in 
the  updated  position. 

The  suggested  algorithm  of  updating  provides  the  validity  of  2  following 
relationships: 


and 


* 

B  =B 


(12) 


w  <  w 

tot  tot 


(13) 
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The  second  step  requires,  in  addition  to  the  current  value  of  r  ,  the  calculation  of 
the  current  mean  curvature  R  of  the  boundary.  Using  current  values  of  T  and  R,  the 
algorithm  requires  the  calculation  of  the  function  %  , 


X  =  2^2 


a 

~R 


(14) 


at  each  point  of  the  current  boundary  3  . 

To  get  the  updated  position  of  the  conductor’s  boundary,  move  each  point  of  the 
surface  in  the  direction  of  the  local  unit  normal  N  .  The  distances  A  of  displacement 
should  be  different  in  different  points  of  the  boundary  (Fig.  3);  it  should  be 
calculated  in  accordance  with  the  relationship 

A  =  -e  (%- x  )  ,  (15) 

where  £  is  a  sufficiently  small  positive  number  and  r  is  the  mean  value  of  the 

J  mean 

of  the  function  %  over  the  surface  3  : 


X 


mean 


area  a  ^ 


(16) 


Fig.  3  Toward  the  2-stage  algorithm  of  computational  analysis  of  the  Rayleigh  instability 

We  call  this  surface  the  upgraded  surface  3*  and  the  domain  B  inside  it  the 

upgraded  domain.  Then,  we  solve  the  electrostatic  problem4,6-9  for  the  upgraded 

* 

domain  and  calculate  the  upgraded  total  energy  W,olal . 

4.  Conclusion 

We  presented  an  algorithm  allowing  one  to  explore  numerically  the  Rayleigh 
problem  dealing  with  destabilization  of  an  electrically  charged  incompressible 
liquid  conductor.  The  algorithm  is  based  on  the  variational  formulation  of  the 
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Rayleigh  problem,  and,  conceptually,  it  allows  one  to  explore  the  under-critical, 
critical,  and  super-critical  regimes  of  charging.  Also,  the  method  permits  further 
extension  for  polarizable  substances. 
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